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Abstract - It has recently been suggested [I] that in a two-level multiplex network, a gradual 
change in the value of the “interlayer” strength p can provoke an abrupt structural transition. The 
critical point p* at which this happens is system-dependent. In this article, we show in a similar 
way as in [2] that this is a consequence of the graph Laplacian formalism used in [T] . We calculate 
the evolution of p* as a function of system size for ER and RR networks. We investigate the 
behavior of structural measures and dynamical processes of a two-level system as a function of p, 
by Monte-Carlo simulations, for simple particle diffusion and for reaction-diffusion systems. We 
find that as p increases there is a smooth transition from two separate networks to a single one. We 
cannot find any abrupt change in static or dynamic behavior of the underlying system. 


Introduction. — In the past several years, single networks have been extensively stud¬ 
ied [3H8] both regarding their structure and also regarding different interactions between 
their nodes. In real world systems, however, there may be more than one type of relation¬ 
ship for the same collection of objects constituting the network. Consider, for example, 
the communication and power networks in a given geographical region. In this case, a fail¬ 
ure in some power station will affect not only the functionality of the power grid, but also 
the routing system for all computers that uses electrical power to sustain its functionality. 
Thus, more recently [9], effort has been applied to the so called “interdependent” or “inter¬ 
connected” networks, meaning a system of two or more networks that are linked together. 
Several articles have been published investigating the properties of these systems mm as 
well as the evolution of dynamical processes on them mm- In their simplest form, they 
consist of two single networks having their nodes connected in a one-to-one configuration 
with each “interlink” having the same strength p (in the range of 0 < p < oo). A first 
question of interest is how does the variation of the parameter p affect the global structural 
properties of the entire system. To this end, one can make use of the Laplacian matrix. 
According to graph theory, given a graph of N nodes, the adjacency matrix A is defined as: 

C = D-A (1) 

where D is a. diagonal N x N matrix with the matrix elements da equal to the vertex 
degree. The lowest eigenvalue of this matrix Ai is equal 0 in the case where the graph is 
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connected. The other eigenvalues contribute to the dynamical processes taking part on such 
a network (e.g diffusion, see m)- The most important role is played by the second smallest 
eigenvalue of the Laplacian, A 2 , which is known as algebraic connectivity |14| . 

Recently, the authors in [1] have used this technique to investigate how the strength 
of the “interlinks” affects the structural properties of a two-level multiplex network. A 
multiplex network is a special case of interdependent networks, where nodes at each layer 
are instances of the same entity. The configuration they used consisted of two coupled 
undirected networks each with the same number of nodes N. Nodes in the fist network were 
connected one-to-one with the nodes of the second network by links of equal strength p. 
Thus, the number of interconnected links is also N. Then, the so called “supra-Laplacian” 
matrix acquires the form: 


^a+pIn -piN A 

V Cb+pIn J ^ ’ 

where /at is the identity matrix of dimension N x N, La and Lb are the Laplacians of 
networks A and B, respectively. The blocks are symmetric matrices of dimension N x N 
(actually this procedure is equivalent to taking the Laplacian matrix of the overall multiplex, 
but now it is in a form which helps to understand the contribution of p). 

The main finding of ref. [T] was that there is a critical value of p, p*, for which the 
“algebraic connectivity” (A 2 ) of the system changes abruptly. For values of p < p*, the 
system behaves as two separate networks, while for p > p* the system behaves like one 
single network. This happens at a value of p for which dX 2 /dp is discontinuous. This 
discontinuity comes as a result of a crossing that occurs between the eigenvalues of the 
supra-Laplacian. In plain terms, above p* the ranking of the eigenvalues changes. Thus, the 
eigenvalue that is initially the second smallest for p < p*, A 2 , becomes larger than another 
eigenvalue (Ai) above p > p*. This is manifested by a discontinuity in the derivative of A 2 
at p*. 

Very recently, Garrahan and Lesanovsky [2] have raised objections at the interpretation 
of the findings in [1] . They pointed out that the eigenvalue crossing is only a consequence of 
the reducibility of the matrix and they applied it for the case of two identical multiplexed 
networks. The characteristic polynomial of the supra-Laplacian matrix L (eq.([5])) is given 
by: 


La + pIn ~ A'^/at 
-piN 


-piN 

Lb + pIn — X^In 


= 0 


li La = Lb, then: 


(3) 


|/:a-a^/jv| |/:a-(a^-2p)j^| = 0 (4) 

From eq.dH), either \La — X^In\ or \La — (A^ — 2p)/Ar| should be equal to 0. For the 
former, we get: Af = Af^ and for the latter Xf = Xf^ + 2p. For a connected system (such 
as network A), 0 = Af*^ < A^*^... < A^'^. Thus, depending on the value of p, we may have 
the following two cases: Af = 0,A2 = Xf^, Af = 2p or Af = OjAf = 2p, Af = Af'^. Thus, 
for a suitable values of X^^ and p, we will have an interchange in the relative position of Af 
and Af and the crossing observed in [1]. 

In the present work, we investigate the effect of the parameter p on various topological 
and dynamical properties of a two-layered multiplex network. Specifically, we want to 
uncover the nature of the p* transition, i.e. is it a sharp, sudden transition or is it a 
smooth one. We calculate numerically the “algebraic connectivity” with respect to interlink 
strength for ER and RR networks of various sizes using standard simulation techniques. We 
investigate the evolution of static structural properties (mean shortest path) and then we 
proceed with various dynamic processes (simple diffusion and reaction-diffusion) to see if 
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we can infer the nature of any structural changes. Furthermore, we investigate the behavior 
of algebraic connectivity as a function of system size. We finish by presenting our main 
conclusions. 

Method. — We used two coupled networks of the same size N and the same statistical 
properties, in a multiplex configuration. The specifications of each configuration are dehned 
in the respective figure captions. In order to specify the dynamical properties of the systems, 
we perform simple particle diffusion for a collection of particles with excluded volume, and 
reaction-diffusion processes of two types of particles {A-\- B model). These are performed 
by standard Monte-Carlo techniques using random walks m- For both processes, we have 
two types of particles, A and B, initially placed in the two different layers, one kind on each 
layer, with a prescribed density p. We let the particles to diffuse. For the simple diffusion, we 
calculate the cumulative number of collisions m between the different types of particles. We 
estimate the value t^oI , which is the time when a certain total number of collisions between 
different types of particles has occurred. For the reaction-diffusion process, where steady 
state condition is imposed, we calculate the mean time required for pN reactions to occur. 
In this model, A reacts with B, but not with another A, and similarly for the B particles. 
We should point out that when two particles react they are removed from the system and 
two new particles (an A and a B) are randomly placed on the layer that corresponds to 
their type. 

We also use the exact enumeration calculation for the case of simple diffusion [16]. In 
the case of a random regular (RR) network, the system is symmetric. Thus, the following 
set of mean-field recursive equations govern the dynamics of the system: 


/Li(t + 1) 
fL2{t + 1 ) 


k +p 

, fL2{t) + 

k +p 


P 

k+p 

P 

k+p 


fL2{t) 

fLl(t) 


(5) 


where fLi{t) is the state of every node on layer LI and fL 2 {t) the state of every node 
on layer L2. In each step t, a node contributes a fraction of to each one of its k 

neighbors on layer i and to the corresponding node on the other layer. We specify 

the time to reach equilibrium (denoted as Tenum) by solving for initial conditions /li(0) = 1 
and /l 2 ( 0 ) = 0, and using the stopping condition ^(/Li(t) - /l2(0) < 10“®- 

The situation is more complicated for the case of ER networks, since they do not have a 
symmetric topology. We consider the mirror case (identical layers), and focus on a node with 
degree k. Denote with ak(t) the containment of this node and with hk{t) the containment of 
its corresponding node (which must have the same degree k). A node with degree k has on 
average ^ ^ nodes with degree k' and from each node a fraction of ak' (t) is passed 

to the node with degree k. The overall contribution should be kJ2k>^i ^ ■ Also, 

a fraction ■^:^bk(t) will be received from the other layer. The same holds for the nodes on 
the other layer. Thus, the following recursive equations hold: 

We solve eq. dH]) for initial conditions afc(O) = 1 and 6^(0) = 0, and using the stopping 
condition \ak{t) — bk{t)\ -i^P{k) < 10“®. Thus, we acquire the time for the system 

to reach the equilibrium, Tenum ■ 
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Furthermore, the change of the structural properties of the system as a function of p 
is investigated using the average shortest path {{Isp)) by varying the interlink strength p 
from nearly zero to 100, using a modified version of the Brandes algorithm m for weighted 
networks. We finally calculate the evolution of p* as a function of N, using the methodology 
proposed in [I]. 

Results and Discussion. — We first examine the evolution of 1/tcoI as a function 
of p, for different cutoff values of m (to being the cumulative number of collisions between 
the different types of particles). In Fig. [1] we present the results for two multiplexed RR 
networks with fc = 4, TV = 1000 and pa = Pb = 0.05, averaged over 100 realizations. We 
observe that the quantity IjTcoi changes smoothly with p. Thus, there is no indication that 
any structural transition can manifest itself through this process. 



Fig. 1: (Color Online), (a) Plot of the inverse of the time, I/tcoI, at which m total collisions between 
A and B particles have been made as a function of p, on RR networks with k — 4, N = 1000 and 
pA = PB ~ 0.05, averaged over 100 realizations. The plot is in log — log scale. Data points are: 
m = 100 (black squares), m = 200 (red circles) and m = 1000 (blue triangles). For all the observed 
cases, the transition is smooth as p increases. 


We next perform the A + B ^ Q reaction-diffusion process and calculate the mean time 
required for all particles to react, Treact- The density of particles A and B is exactly the 
same (p = 0.25) and remains hxed throughout the process (steady state). All A particles 
are initially placed on one layer and all B particles on the other. We calculate the time 
Treact which is the time it takes to have pN reactions. In Fig. [21 we plot the first derivative 
of \/Treact as a function of p. As p increases, this quantity tends to 0, meaning that Treact 
is independent of p for large p values. If the system were to undergo an abrupt structural 
transition for a specific value of p, we should have observed qualitatively this sharp change 
in the vicinity of such a point. 

Finally, we use the exact enumeration method, for the case of the diffusion process. In 
Fig. m we plot the inverse of the equilibration time (l/Tenum) as a function of p for a 
setting of two identical multiplexed RR networks (/c = 4, TV = 1000) and ER networks 
LFigTO 
in d]; 


((fc) = 5,TV = 1000). In 
with the method used 


a peak appears at p ~ 4, which is not the value calculated 
which was p* ~ 0.41. The same behavior is recovered when 


fc = 8 for RR networks (data not shown). The results are accompanied by the numerical 
simulation of the process, and as we see they are in excellent agreement. A qualitative 
explanation can be given here: when p < fc, it happens that p/fc + p < k/k + p, and the 
majority of particles remain for longer times in the starting layer. As p increases particles 
tend to move to the other layer, thus promoting the decrease of the equilibration time. 
When p > fc, we have that p/fc + p > fc/fc + p, and thus the inverse procedure takes place. 
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Fig. 2: Plot of the derivative of the inverse characteristic time (1/rreact) as a function of p, for two 
identical ER networks with (fc) = 5 and N — 1000. 


The case is similar for coupled ER networks, see Fig. |3(b)[ except for the fact that the peak 
is not sharp since in this case there is an inherent asymmetry in the degree distribution. 

These results indicate that there is no manifestation of any sudden structural change on 
common dynamical processes on multiplex networks. We further proceed by investigating 
network properties of the multiplex systems. We search to see if there is some structural 
modification in the topological properties of the network by investigating the “average short¬ 
est path” {{Isp})- We performed calculations on a random-regular multiplex network with 
k = 8 and N = 1000, using the algorithm proposed in m for the case of weighted graphs. 
Regarding the weights, we used a value of 1 for all intralayer connections and p for the 
interlayer ones. We define the distance between two nodes connected by an interlink as the 
inverse of its strength p. The results for (Isp) are shown in Fig. |4l where we observe that 
as p —>■ 0, (Isp) oo because 1/p oo. Furthermore, when p <C 1, (Isp) is controlled by 
the maximum shortest path between a source node in network A(B) and a target node in 
network B(A). As p increases, 1/p decreases and so does (Isp)- For p 2> 1, the multiplex be¬ 
haves as a single network. We observe that, by varying the value of p, this quantity changes 
in a continuous way. Thus, no abrupt transition can be seen here from the data as well. 

We now apply the methodology proposed in [1] (see also eq. ([2])) to investigate the 
behavior of the algebraic connectivity, A 2 , as a function of the size of the layers. In Fig. [Sj 
we plot the evolution of p* as a function of system size N for two types of networks: ER 
networks with (k) = 5 (Fig 5(a)) and random regular with k = 8 (Fig 5(b)). We observe a 
clear dependence of p* on the size N of the layers. Also, there is a fast decay on the value 
of p*, but we cannot determine its nature (i.e. being power-law or exponential), due to the 
small range of N. 


Conclusions. — In this paper we have tested structural as well as dynamic properties of 
multiplex networks. We used various types of networks such as coupled RR and ER (identical 
or not) networks and various different methodologies. We have observed a continuous change 
of the structural properties as a function of p, as shown by the behavior of the average 
shortest path {{Isp))- We have performed diffusion and reaction-diffusion processes on the 
above multiplex configurations studying suitable characteristic times associated with the 
structural properties of the system. In all cases we verified that changes take place in a 
continuous manner. Also, the point where such changes occur does not coincide with the p* 
calculated with the method proposed in [1]. 

Ackno-wledgments. — Results presented in this work have been produced using the 
European Grid Infrastructure (EGI) through the National Grid Infrastructures NGI-GRNET^ 
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(a) 



(b) 

Fig. 3: (Color Online). Plot of the inverse of the equilibration time {1 jTe.nii.m) as a function of p, 
for the case of (a) two coupled random-regular networks with fc = 4 and N = 1000 and (b) two 
coupled identical ER networks with {k) = 5 and N = 1000. In (a), there is a peak at p ~ 4. In 
both figures, there is a change in the behavior for a value of p. However the point at which this 
occurs does not coincide with the p* calculated with the method proposed in [I] . 
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Fig. 4: Plot of the “average shortest path length” {{lap}) as a function of p for a multiplex of two 
random regular networks with k = 8 and N = 1000 each. It is obvious that there is no abrupt 
transition of this quantity. 
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Fig. 5: Plot of the value of the transition point of the algebraic connectivity, p*, as a function of 
the size of each layer N, on (a) two coupled ER networks with {k) = 5 and (b) two coupled RR 
networks with fc = 8. Data for both cases are averaged over 100 rnns for large systems and up to 
5000 for small ones. It is evident that p* is system and size dependent. 
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